function FVlinearadvectionFOU2D 
% File name: FVlinearadvectionFOU2D.m 
% Description:  Solves the 2D PDE, dU/dt + div(H) = 0 
% using the FOU finite volume scheme. 
% This equation corresponds to the finite volume form: 
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0. 
% 
% In this case the equation is the linear advection equation so: 
% 
% flux density: H = vU  
% velocity: v = vx i + vy j = <vx,vy> where vx and vy are constant. 
% 
% The program is written in a structured way using subfunctions so that it 
% can easily be modified for other equations, schemes and meshes. 
% subfunctions: freadmesh, fcellarea, fsidevectors, finitialu, 
%               fIflux, fghostcells, fcalcdt, fplotresults, fdisplay, 
%               fdrawmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Verification:  (on uniform Cartesian mesh) 
% for NI=61,NJ=81,dx=2,dy=4, initial profile centre is at (61,82)  
% contour plots of numerical agree with exact for:  
% vx=10,vy=0, t=1.04, new centre (71, 82) is correct 
% vx=0, vy=10, t=1.12, new centre (61, 94) is correct 
% vx=10,vy=10, t=1.04, new centre (71, 94) is correct 
%  
% Comments:  Although program appears to work the numerical scheme is highly 
%            dissipative. 
  
%% 
clc; clf; clear all; 
runtime=1.0;     % target runtime [s] 
t=0.0;           % time [s] 
tlevel=0;        % time level 
vx=10;   % velocity component in x direction [m/s] 
vy=10;   % velocity component in y direction [m/s] 
%vx=0;vy=10; 
v=[vx,vy];       % velocity vector <vx,vy> 
[x,y,xcen,ycen,solid]=freadmesh; % gets coordinates of the 2D computational mesh  
                     % cell centres and solid flags (excluding ghost cells) 
disp('read in mesh') 
[m,n]=size(x);  
NI=m-1; % number of computational cells in i direction. 
NJ=n-1; % number of computational cells in j direction. 
A=fcellarea(x,y);    % compute and store cell areas in an NIxNJ array 
disp('calculated cell areas') 
S=fsidevectors(x,y); % compute and store cell side vectors in an 
                     % (NI+2)x(NJ+2)x4x2 array which contains ghost cells. 
disp('calculated cell side vectors') 
u=finitialu(xcen,ycen,0,v); % put initial cell centre values of u in NIxNJ array 
uinitial=u;          % store for plotting initial conditions 
disp('inserted initial conditions') 
% 
% Append extra cells around arrays to store any ghost values. 
u=finsertghostcells(u); % u is now (NI+2)x(NJ+2) 
unext=zeros(size(u));   % (NI+2)x(NJ+2) array for u values at next time level 
A=finsertghostcells(A); % A is now (NI+2)x(NJ+2) 
%% 
disp('time marching starts')    
while(t<runtime)  
 u=fbcs(u);  % Implement boundary conditions using ghost cells.  
             % u is (NI+2)x(NJ+2) and each computational cell is  
             % indexed i=2 to NI+1, j=2 to NJ+1. 
 dt=fcalcdt(A,S,v); % Finds stable time step for each iteration. 
 t=t+dt  % update time  
 tlevel=tlevel+1;   
 for j=2:NJ+1  
   for i=2:NI+1  
       IH=fIflux(v,u,i,j); % gets interface fluxes for cell (i,j) 
       IH1=[IH(1,1),IH(1,2)]; % interface flux vector for side 1 
       IH2=[IH(2,1),IH(2,2)]; % interface flux vector for side 2 
       IH3=[IH(3,1),IH(3,2)]; % interface flux vector for side 3 
       IH4=[IH(4,1),IH(4,2)]; % interface flux vector for side 4 
       % 
       s1=[S(i,j,1,1),S(i,j,1,2)]; % side vector for side 1 
       s2=[S(i,j,2,1),S(i,j,2,2)]; % side vector for side 2 
       s3=[S(i,j,3,1),S(i,j,3,2)]; % side vector for side 3 
       s4=[S(i,j,4,1),S(i,j,4,2)]; % side vector for side 4 
       % FV scheme 
       % compute total flux out of cell (i,j) 
       totalfluxout=dot(IH1,s1)+dot(IH2,s2)+dot(IH3,s3)+dot(IH4,s4); 
       unext(i,j)=u(i,j)-(dt/A(i,j))*totalfluxout;  
   end % of i loop 
 end % of j loop         
 u(2:NI+1,2:NJ+1)=unext(2:NI+1,2:NJ+1);  % update u values 
 % 
end  % of while loop 
disp('time marching ends')  
uexact=finitialu(xcen,ycen,t,v); % exact solution at time t 
u=u(2:NI+1,2:NJ+1); % extract final computational values 
fdrawmesh(x,y,solid) % draws mesh if fplotresults is commented out 
%fdisplay(u); % print results as a matrix for a small mesh if needed 
fplotresults(xcen,ycen,uinitial,u,uexact) % plot results: comment in or out 
disp('finished main program, see figures') 
end % FVlinearadvectionFOU2D 
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [x,y,xcen,ycen,solid]=freadmesh 
% The mesh is structured and has NI cells in the i direction and  
% NJ cells in the j direction (so each cell has 4 sides). 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively which are both (NI+1)x(NJ+1).  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  xcen and ycen are NI by NJ arrays 
% which hold cell centres. solid is an NI by NJ array which  
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise 
% solid(i,j)=0. 
% 
% % test using diamond 
% x(1,1)=0; y(1,1)=-1; 
% x(2,1)=1; y(2,1)=0; 
% x(2,2)=0; y(2,2)=1; 
% x(1,2)=-1; y(1,2)=0; 
% 
% Test using regular Cartesian mesh starting at (0,0).  i direction is 
% x direction. 
% 
% NI=61,NJ=41,dx=2,dy=4 makes mesh centre (61,82) 
NI=61; 
NJ=41; 
dx=2; 
dy=4; 
x=zeros(NI+1,NJ+1); % allocate correct sized array 
y=zeros(NI+1,NJ+1); % allocate correct sized array 
xcen=zeros(NI,NJ);  % allocate correct sized array 
ycen=zeros(NI,NJ);  % allocate correct sized array 
solid=zeros(NI,NJ); % allocate correct sized array 
% create cell vertices 
for i=1:NI+1 
    for j=1:NJ+1 
        x(i,j)=(i-1)*dx; 
        y(i,j)=(j-1)*dy; 
    end % of j loop 
end % of i loop 
nonCartesian=1; % 1 for non-Cartesian mesh 
if (nonCartesian==1) 
   disp('mesh is non-Cartesian') 
   % Add 'random' numbers to mesh coordinates to produce a non-Cartesian mesh 
   % ensure that the abs of these numbers are less than min(dx/2,dy/2). 
   % Create non-Cartesian mesh 
   for i=1:NI+1 
    for j=1:NJ+1 
        x(i,j)=x(i,j)+(dx/6)*cos(1.727*(i+j)); 
        y(i,j)=y(i,j)+(dy/5)*sin(2.46723*i*j); 
    end % of j loop 
   end % of i loop 
else 
   disp('mesh is Cartesian')  
end % end of if statement create mesh 
% 
% find cell centres by averaging coordinates of vertices 
for i=1:NI 
    for j=1:NJ 
        xcen(i,j)=(x(i,j)+x(i+1,j)+x(i+1,j+1)+x(i,j+1))/4; 
        ycen(i,j)=(y(i,j)+y(i+1,j)+y(i+1,j+1)+y(i,j+1))/4; 
    end % of j loop 
end % of i loop 
end % freadmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function A=fcellarea(x,y)   
% Verification:  correct for single diamond 
% cell area = half the modulus of cross product of cell diagonal vectors. 
% These vectors are 2D so append 0 k component to make 3D for cross product 
[m,n]=size(x); 
NI=m-1; NJ=n-1; 
A=zeros(NI,NJ); 
for i=1:NI 
    for j=1:NJ 
        d1=[x(i+1,j+1)-x(i,j),y(i+1,j+1)-y(i,j),0]; 
        d2=[x(i,j+1)-x(i+1,j),y(i,j+1)-y(i+1,j),0]; 
        d=cross(d1,d2);    
        A(i,j)=0.5*sqrt(dot(d,d)); 
    end % of j loop 
end % of i loop 
end % fcellarea 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function S=fsidevectors(x,y) 
% Verification: correct for single diamond shaped cell. 
% 
% Finds the 4 side (S) vectors for each cell.   
% Cell sides are numbered 1 to 4 anti-clockwise.  For cell (i,j) side 1  
% is the vector s1 from (x(i,j), y(i,j)) to (x(i+1,j),y(i+1,j)). 
% if s1=<a,b> then the corresponding side vector S1=<b,-a>. 
% The side vector array has ghost cells around it so that side vector  
% indices (i+1,j+1) refer to cell (i,j). 
% 
% Note that vectors are 1 by 2 arrays so care must be taken because side 
% vectors are stored in higher dimensional arrays. 
% 
[m,n]=size(x); 
NI=m-1; NJ=n-1; 
S=zeros(NI+2,NJ+2,4,2); 
for i=1:NI 
    for j=1:NJ 
        s1=[x(i+1,j)-x(i,j),y(i+1,j)-y(i,j)]; 
        S(i+1,j+1,1,1)=s1(2); 
        S(i+1,j+1,1,2)=-s1(1); 
        % 
        s2=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)]; 
        S(i+1,j+1,2,1)=s2(2); 
        S(i+1,j+1,2,2)=-s2(1); 
        % 
        s3=[x(i,j+1)-x(i+1,j+1),y(i,j+1)-y(i+1,j+1)]; 
        S(i+1,j+1,3,1)=s3(2); 
        S(i+1,j+1,3,2)=-s3(1); 
        % 
        s4=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)]; 
        S(i+1,j+1,4,1)=s4(2); 
        S(i+1,j+1,4,2)=-s4(1);     
    end % of j loop 
end % of i loop 
% [S(2,2,1,1),S(2,2,1,2)] 
% [S(2,2,2,1),S(2,2,2,2)] they are correct for a diamond 
% [S(2,2,3,1),S(2,2,3,2)] 
% [S(2,2,4,1),S(2,2,4,2)] 
end % fsidevectors 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function u=finitialu(xcen,ycen,t,v)    
% Inserts exact u values (and initial conditions for t=0)  
[NI,NJ]=size(xcen); 
u=zeros(NI,NJ);  % create correct sized array (without ghost values) 
vx=v(1); 
vy=v(2); 
% 
% Initial 2D Gaussian function based on the centre (xc,yc) of the  
% computational domain. 
xmax=max(max(xcen));  % max x value in mesh 
ymax=max(max(ycen));  % max y value in mesh 
xmin=min(min(xcen));  % min x value in mesh 
ymin=min(min(ycen));  % min y value in mesh 
xc=(xmax+xmin)/2;  % approx x coord of mesh centre 
yc=(ymax+ymin)/2;  % approx y coord of mesh centre 
% 
cutoff=min(xmax-xc,ymax-yc)/2;  % cut off radius for Gaussian 
for i=1:NI 
    for j=1:NJ 
        x=xcen(i,j)-vx*t;  
        y=ycen(i,j)-vy*t; 
        d=sqrt((x-xc)^2+(y-yc)^2); % distance from centre 
        if (d<cutoff) 
           u(i,j)=exp(-0.01*d^2); 
        else 
           u(i,j)=0; 
        end 
    end % of j loop 
end % of i loop 
end % finitialu 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function IH=fIflux(v,u,i,j) 
% Calculates the 4 interface fluxes on the 4 sides of cell (i,j).  
% Cell sides are numbered 1 to 4 anti-clockwise.  For cell (i,j) side 1  
% is the vector s1 from (x(i,j), y(i,j)) and to (x(i+1,j),y(i+1,j)). 
% Interface fluxes are denoted IH1, IH2, IH3 and IH4 and are 2D vectors. 
% Flux H depends on U, i.e. H = H(U).  
% For the 2D linear advection equation, H(U) = <vx U, vy U>. 
%  
% There is considerable choice for interface flux estimation: different  
% choices give different FV schemes. 
%  
%% The following scheme simply takes: 
%  IH1=H(u(i,j-1)) = v u(i,j-1)  
%  IH2=H(u(i,j)) = v u(i,j) 
%  IH3=H(u(i,j)) = v u(i,j) 
%  IH4=H(u(i-1,j)) = v u(i-1,j) 
% This is an upwind scheme for vx > 0 and vy >0 and corresponds to the  
% FOU finite difference scheme on a Cartesian mesh. 
% 
IH=zeros(4,2); % array to store all 4 interface flux vectors for a cell 
% 
vx=v(1); 
vy=v(2); 
% 
IH(1,1)=vx*u(i,j-1);  % component of IH1 in the x direction 
IH(1,2)=vy*u(i,j-1);  % component of IH1 in the y direction  
IH(2,1)=vx*u(i,j);    % component of IH2 in the x direction 
IH(2,2)=vy*u(i,j);    % component of IH2 in the y direction 
IH(3,1)=vx*u(i,j);    % component of IH3 in the x direction 
IH(3,2)=vy*u(i,j);    % component of IH3 in the y direction  
IH(4,1)=vx*u(i-1,j);  % component of IH4 in the x direction 
IH(4,2)=vy*u(i-1,j);  % component of IH4 in the y direction 
end % fIflux 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function outarray=finsertghostcells(inarray) 
% Assuming that ghost cells are needed around the whole domain an mxn 
% array is embedded into (m+2)x(n+2) array of zeros.  Hence computational  
% indices go from i=2 to m+1 and j=2 to n+1. 
[m,n]=size(inarray); 
dummy=zeros(m+2,n+2); 
dummy(2:m+1,2:n+1)=inarray; 
outarray=dummy; 
end % finsertghostcells 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% function outarray=fremoveghostcells(inarray) 
% % Removes ghost cells so that a mxn array becomes (m-2)x(n-2) 
% [m,n]=size(inarray); 
% outarray=inarray(2:m-1,2:n-1); 
% end % fremoveghostcells 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function u=fbcs(u) 
% Implements boundary conditions using ghost cells index by i=1, i=NI+2 
% j=1, j=NJ+2. 
[m,n]=size(u); 
NI=m-2; NJ=n-2; 
% 
% Dirichlet: u=0 everywhere. 
% Don't need to do anything as u values in ghost cells are already zero. 
end % bcs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function dt=fcalcdt(A,S,v) 
% Verification: 
% Finds allowable time step dt using heuristic formula (which makes sense 
% 0<F<1 is a tunable safety factor. 
[m,n]=size(A);  % A and S now enlarged to include ghost cells 
NI=m-2; NJ=n-2;  
dt=9999999;  % large value of dt to start 
F=0.4; 
 for i=2:NI+1 
    for j=2:NJ+1 
        % put side vectors into 1x2 arrays   
        s2=[S(i,j,2,1),S(i,j,2,2)]; 
        s3=[S(i,j,3,1),S(i,j,3,2)]; 
        % 
        v2=abs(dot(v,s2)); % side 2 
        v3=abs(dot(v,s3)); % side 3 
        dum=A(i,j)/max(v2,v3); 
        % take smallest value of dt 
        if (dum<dt) 
           dt=dum; 
        end 
    end % j loop 
 end % i loop 
dt=F*dt; 
end  % fcalcdt 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fdrawmesh(x,y,solid) 
% Description:  Draws a structured 2D finite volume mesh and fills in any 
% solid cells.  
% Date structures: 
% The mesh has NI cells in the i direction and NJ cells in the j direction. 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively which are both NI+1 by NJ+1.  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  solid is an NI by NJ array which  
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise 
% solid(i,j)=0. 
%  
[m,n]=size(x); 
NI=m-1; % number of cells in i direction 
NJ=n-1; % number of cells in j direction 
% 
% Plot the mesh 
hold on   % keeps all plots on the same axes 
% draw lines in the i direction 
for i=1:NI+1 
    plot(x(i,:),y(i,:)) 
end 
% draw lines in the j direction 
for j=1:NJ+1 
    plot(x(:,j),y(:,j))    
end 
title('computational mesh') 
xlabel('x') 
ylabel('y') 
% Fill in solid cells 
for i=1:NI 
    for j=1:NJ 
        if (solid(i,j)==1) 
            solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)]; 
            solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)]; 
            fill(solidx,solidy,'k')       
        end 
    end % of j loop 
end % of i loop 
end % fdrawmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fplotresults(x,y,uinitial,u,uexact) 
% Plots results on a structured mesh - even works for non-Cartesian! 
% 
subplot(3,1,1), contour(x,y,uinitial) 
title('Contour plot of initial conditions') 
xlabel('x [m]') 
ylabel('y [m]') 
% 
subplot(3,1,2), contour(x,y,u) 
title('Contour plot of numerical solution') 
xlabel('x [m]') 
ylabel('y [m]') 
% 
subplot(3,1,3), contour(x,y,uexact) 
title('Contour plot of exact solution') 
xlabel('x [m]') 
ylabel('y [m]') 
% 
end % fplotresults 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fdisplay(in) 
 % displays matrices 'correctly' - WORKS 
 % Using index i for x which are columns and j for y which are rows in  
 % reverse order so this routine displays them in this fashion 
 [NX,NY]=size(in); 
 temp=transpose(in); 
 out=temp(NY:-1:1,:); 
 disp('printing as a matrix') 
 out 
end % fdisplay 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 

